Macroscopic Coulomb blockade in large Josephson junction arrays 
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We investigate theoretically transport properties of one- and two-dimensional regular Josephson 
junction arrays (JJAs) in an insulating state. We derive the low-temperature current-voltage char- 
acteristics (the I-V dependencies) for the current mediated by the Cooper pair transfer across the 
system. In the case where the screening length Ac associated with the capacitance of the islands 
to the ground is much larger than the island's size d, we find that transport is governed by the 
macroscopic Coulomb blockade effect with the gap Ac well exceeding a single island charging energy 
Ec- In the limit of Ac ^ L, where L is the linear size of the array, the gap establishes the dependence 
on the array size, namely, Ac — Ec{L/d) in ID and Ac — -Ecln(L/d) in 2D arrays. We find two 
transport regimes: at moderate temperatures, Ec < ksT < Ac, the low bias transport is thermally 
activated with the resistance R oc exp{To/T) where the activation energy ksTo = Ac. At ultra-low 
temperatures, ksT < Ec, a JJA falls into a superinsulating state with R oc exp[{Ac/ Ec) exp(i?c/r)]. 

PACS numbers; 05.60.Gg,74.81.Fa,73.63.-b 



I. INTRODUCTION 

A theoretical and experimental study of large regular 
Josephson junction arrays (JJAs), the systems comprised 
of small superconducting islands connected by Joseph- 
son junctions, has a long historji^i^i^i^i^iii^i^iiSiiiii^ii^, 
ggQi4ji5 £qj. review. The remarkable feature of these 
systems is that they experience a superconductor-to- 
insulator transition (SIT) as the Coulomb charging 
energy of a single island, Ec (i.e. an energy cost 
to place a Cooper pair on such an island) compares 
to the Josephson coupling energy, Ej measuring the 
strength of the phase coupling between the supercon- 
ducting islands comprising a JJA In arrays 
near the critical balance between Ec and Ej, a SIT 
can be induced by the magnetic field^^. Studies on 
the superconductor-insulator transition in thin granular 
films, which are often modeled as JJAs, revealed the 
similar behavioria^iSi^iSS^iSS^^^^^^^. Even 
more remarkably, critically disordered homogeneous su- 
perconducting films exhibited all the wealth of phe- 
nomena related to superconductor-insulator transi- 
tion characteristics to granular superconducting sys- 
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brought about the idea that strong disorder induces 
an inhomogeneous spatial structure of isolated super- 
conducting islands in thin homogeneously disordered 
films^li^Si^i^iSi^. Numerical simulations of the ho- 
mogeneously disordered superconducting films confirmed 
that indeed in the high-disorder regime, the system 
breaks up into superconducting islands separated by an 
insulating se a^^i^^'^"^ . Recent scanning tunneling micro- 
scope measurements of the local density of states in thin 
TiN films^* offered strong support to this hypothesis. 

All the above together shows that the Josephson junc- 
tion arrays offer a useful generic model that captures 
most essential features of the superconductor-insulator 



transition in a wide class of systems ranging from ar- 
tificially manufactured Josephson junction arrays to su- 
perconducting granular systems and even the homoge- 
neously disordered superconducting films (see also re- 
view^^) and allows for consideration of all of them on 
the common ground. Indeed, recent theoretical results 
describing insulating behavior of regular JJAs, appeared 
to be in a striking quantitative accordance with the ex- 
perimental findings in TiN and InO superconducting dis- 
ordered films^i^. 

The common "working tool" for experimental study 
of the SIT is the measurements of transport charac- 
teristics of the systems in question. Altering vari- 
ous parameters of the system, such as tunnel resis- 
tance and transmittance in JJAs, conditions of deposi- 
tion, chemical composition, and thickness of the films, 
and applied magnetic field, one can drive the system 
directly from the superconducting to insulating state. 
The transition is observed as a set of the fan-shaped 
temperature d epend ences of the resistance R(T), see 
Refs. jr 3"4,6',9' 12 13 19 21 24 27 28 29,30','3lll3i|33l3il3i, 
^^mMMLMiMMd5Mi£LMMi5Q,l with the acti- 
vation behavior 



R{T) cx e^°/^ 



(1) 



on the insulating side of the transition, see Refs. [sIlTllslfl^. 
!T3J9,35,36,37,38^43^47,48^49.5(|, where Tq is the activa- 
tion temperature. Most pronounced features of this tran- 
sition manifest themselves in the current- voltage charac- 
teristics (the I-V curves). On the superconducting side 
a system has very low resistance at low currents followed 
by jump in resistivity when current exceeds the critical 
value. On the insulating side the I-V characteristics show 
a mirror behavior: extremely high resistance at low volt- 
ages and abrupt jump in the conductivity at the thresh- 
old voltage 14-, see Refs. [iliiiliQiiiHiSiiiiSii, 

HOHmillil . Yet the most startling observations that 



come from the insulating side of the transition are: (i) 
the size dependent activation energy£9; (ii) hyperactiva- 
tion temperature behavior of resistivity at ultra- low tem- 
peratures in JJAs^'^; (iii) size dependent threshold volt- 
ageii^i^; (iv) peculiar interrelated magnetic field depen- 
dencies of the activation energy and threshold voltage in 
JJAsfi and superconducting films^iiS. 

The above striking findings called for a theory capable 
of quantitative description of the accumulated wealth of 
the experimental results within a unified picture. The 
preceding publications^"-^^ offered such a description in 
a framework of the Cooper pair transport in large JJAs in 
the insulating region, where Ec ^ Ej, Ej = fi/c/2e {Ic is 
the critical current of a single Josephson junction). This 
paper is the extended version of earlier publications^!^, 
presenting the details of the derivation of the current- 
voltage characteristics. 

We briefly summarize our results: The insulating be- 
havior of the large Josephson junction array is governed 
by the macroscopic Coulomb blockade effect with the 
Coulomb blockade activation energy 

_ ( E,[A/i2d)] , li? array 
~ \ Ec\n{A/d) , 21? array ' 

where A — min{L, Xc}, L is the size of an array, d 
is the size of the elemental unit of JJA, and Ac is 
the screening length related to the capacitance to the 
ground. Importantly, this activation energy can be 
much larger than the single junction charging energy 
Ec- In the two-dimensional array the charge binding- 
unbinding Berczinskii-Kosterlitz-Thouless like transition 
takes place at T = Tgi ~ Ec/ks separating the insu- 
lating phase existing in the interval Ec < ksT < Ac 
and exhibiting the thermally activated resistivity (1) with 
ksTo = Ac, and the superinsulating state at T < Ec/ks, 
with R oc exp[(Ac/£'c) exp(i?c/fcB2^)]- In the one dimen- 
sional arrays we expect a crossover between these two 
states. 

The paper is organized as follows: in Section II we 
introduce the model and present the general equation for 
the dc I-V dependence, expressing it through the time- 
dependent correlation function K{t) of superconducting 
order parameter phases of the whole system. In Sections 
III and IV we discuss the properties of K{t) and the 
corresponding I-V characteristics above and below Tsi 
respectively. Section V presents the discussion of the 
obtained results. 



II. MODEL AND CVC CALCULATION 

The current-voltage characteristics of Josephson sys- 
tems in an insulating state have been a subject of intense 
discussions over the decades. Single junctions were con- 
sidered and the behavior of two-junction 
systems was studied i n^^i^°i^^ . Only a few works were 
addressing the Cooper-pair current in large JJA a^°i^^'^^ . 
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FIG. 1: Sketch of the considered array geometries. An exter- 
nal current I is injected from the left through the electrode 
having the superconducting phase Xl and extracted through 
the right electrode with the phase Xh- Upper panel: One- 
dimensional array of N superconducting islands (squares) 
connected by two Josephson junctions (crosses) to neigh- 
bors corresponding to experimental system of Ref. 5] . Lower 
panel: Two-dimensional A'^ x M Josephson junction array. 



In this paper we will discuss Cooper pair transport 
(Josephson current) and the corresponding I-V char- 
acteristics in large JJAs, leaving the calculation of the 
quasiparticle contribution to the forthcoming publica- 
tion. This implies, in particular, that we neglect interac- 
tions of the internal phases with the thermal bath, since 
the latter is equivalent to switching on quasiparticle cur- 
rent. Let us consider N x M superconducting islands 
comprising a one- (Af = 1) or two-dimensional array 
closed by a small (as compared to the quantum resis- 
tance for Cooper pairs Rcp — h/Ae^ ~ 6.45 kfi) external 
resistance, Rext^ see Fig. I. Note, that in this kind of 
circuits with a small external load resistance both, a sin- 
gle Josephson junctionS*'^^'^^ and two Josephson junc- 
tions in series— are always in a superconducting state 
irrespectively to relation between Ec and Ej . We are in- 
terested in a low-temperature transport, T <^ Tc, where 
Tc is the critical temperature of a single superconduct- 
ing island, and, therefore, we can neglect the fluctua- 
tions of the amplitude of the order parameter. We as- 
sign the fluctuating order parameter phase Xij (t) to the 
{i, j}-th superconducting island (see Fig. I). Josephson 
relation connects these phases with the fluctuating volt- 
age drops between the adjacent islands. We denote the 
phases of the left- and right leads as Xt(0 ^-nd Xn{t), cor- 
respondingly. The finite voltage V applied to a JJA gen- 
erates the alternating Josephson currents proportional to 
Ej sin{Vt -\- {Xij}{t))- The dc component of the Joseph- 
son current results from the time averaging of the ac 
currents and is thus determined by the correlations of 
the time-dependent fluctuations of the Josephson phases 
{Xij}{t) across the array. 

We discriminate between the Josephson coupling ener- 
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gies of intrinsic junctions Ej and the Josephson coupling 
energies Ej between leads and first and A^-th rows. Sin- 



J 



gling out the terms containing the leads phases Xl 
Xr explicitly we present the array Hamiltonian as: 



and 



M 



H = Ho + Hbath + Hm{Xl ~ Xh} + -T^Y.^{Xlj{t)-XL f + {XNj{t)-Xnf] 

M . 

- Ej^l cos [xL{t) - xij(0] + cos [xiiit) - XNjit)] 



Here 



-ffn = I -T^iXi] - Xki - 2eVij-ki/hy - Ej cos{xtj - Xki) 



E 



(3) 



(4) 



r 



the brackets (ij, kl) denote summation over the pairs of 
adjacent junctions, and the last term in Q represents 
the self-charge energies of superconducting islands. We 
introduced here the dc voltage drops on the junctions 
Vij^ki- The charging energies Ec and Eco are determined 
by the junction capacitance C and capacitance to the 
ground Co, as Ec = 2e^/C and Eco — 2e^/Co, respec- 
tively. The Hhath is the Hamiltonian characterizing the 
thermal bath, which can be modeled as a set of harmonic 
oscillators with coordinates ^ir^^, 



H, 



bath 



E 



2A/,: 



I.e. 



(5) 



where uji and Mi are the frequency and mass of harmonic 
oscillators. The Hint term in ([3]) describes a bilinear 
coupling of phases on the leads to the thermal heat bath 
as 



(6) 



where Ai are the coupling constants. In the cir- 
cuits presented in Fig. 1, the coupling constants 



J 



Ai are determined by the external resistance^i^i 
J2i ^1 / i-^i^i) — Rcp/Rext- Since we consider Joseph- 
son currents only, we will neglect the phase coupling to 
the thermal heat bath in the internal part of the array: 
dissipation on the internal islands implies the presence 
of the current of quasiparticles. This means that we 
can treat the evolution of the internal phases as a non- 
dissipative quantum dynamics. 

All the above implies that the dynamics of the phases 
in the leads Xr ^-nd Xl interacting with the heat bath 
differs from that of the internal phases, i.e. the phases 
on superconducting islands Xij- Namely, whereas the 
phases in the leads are to be treated as classical dy- 
namic variables satisfying the Langevin stochastic equa- 
tion, the phases in the internal superconducting islands 
are the quantum-mechanical variables with the dynam- 
ics described by the quantum-mechanical Hamiltonian 
Hq. Moreover, the phases in the leads and the intrin- 
sic phases interact through the Josephson coupling terms 
[the last two terms in the array Hamiltonian ([3])]. Now, 
shifting all the phases over (xt + Xr)/2, i-e replacing 
Xi3 Xij - {Xl + Xfl)/2, we obtain: 



^2 ^2 J" 

H = Ho + Hbath + HintixL - Xr} + E(^« ^ Xl - Xl] + XNjf + 'Y^iXij - XNjf 



2 \Ec Eco 



i=i 



1 1 \ . 

- + — - 2Ej ^ cos 0, cos 



Xr - Xl- Xl]{t) +XNjit) 



(7) 



r 



where a new mdependeni variable — (xij+XAfj)/2 had array is calculated as /s = {dH/d[xR — XL])- Considering 
been introduced. The dc Josephson current through the Xr ~ Xl as a parameter, one can apply the Hellmann- 
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Feynman theorem2§. and reduce this to^SiZi 



The term in the total Hamiltonian depending exphcitly 



Is{V) 



d[xR - Xl 



on the variables is 



(8) 



2 \E. 



E, 



cO 



M 

E 

i=i 



M 

2Ej ^ cos X cos 



r 



Xr - Xl- Xljit) + XNj{t) 



(9) 



As a first step, we perform averaging over In the 
insulator regime, which we address here, Ec,Ecq ^ Ej, 
and the averaging procedure is carried out by making use 
of the perturbation theory with respect to Ej/E^ [the 
last term in Eq. 0]. Such a procedure is similar to the 
one used in the case of the Cooper pair two-junctions 
transistor^^'''-. In the zero order perturbation theory 



{cos = 0. Indeed, {cos cx {n\ cos 4'j\n) = 0, 

where \n) = Y\j^^{1/\/2tt) exp{i(j)jn) are the wave func- 
tions corresponding to the first term in the Hamiltonian 
©, i.e. the wave functions in a zero order of the pertur- 
bation theory. The wave functions in the first order of 
the perturbation theory in Ej/E^ become^^ 



\iJ„^\n)- 2Ej cos 



Xr - Xl -Xi]{t) +XNj{t) 



M 



EE 

j = l m 



p(0) p(0) 
-tin — -C/m 



■|m) , 



(10) 



where 



EcEca'n? 



2{Ec + Eco) 

are the energy levels obtained in a zero order of the perturbation theory Therefore, in the first order one finds 



cos 



Xr-Xl -Xij{t) + XNj{t) 



COS( 



Ej {Ec + Ec 



PT('t>.) E,E,o ^ (n2-l/4) 

n— — oo ^ ' ' 



E 



Xr -Xl- Xi]{t) +XNj{t) 



(11) 



Here, we used that for finite temperatures the quantum-mechanical distribution of 0^ is determined by the equilibrium 
density matrix corresponding to a first term in the Hamiltonian ([9]), i.e. 

pA^j) ^ ^e-^|n)(n|. 

n 

Further, to deal with less awkward formulas, we consider the case C 3> Co (the general situation with the arbitrary 
relation between C and Co is straightforwardly recovered as needed) and introduce the parameter 



E 



n— — oo 



(2^2-1/2) ■ 

The effective Hamiltonian H^f / depending on the variables Xr — Xl and Xij assumes the form: 

^eff — Efo + Htath + HintiXR ~ Xl} + El* , 

where 



M 



= E -Xl- Xlj{t) + XNj{t)] , 



E. 



(12) 



(13) 
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and, accordingly, the dc component of the non-dissipative Cooper pair Josephson current across the array is 

dH* 



Is{V) 



IcM hm - / dt{ sin [xr - Xl 

he r^oc T Jq \ 



Xlj{t) +XNjit)] 



(14) 



where, the brackets (...) stand for a quantum mechanical averaging over phases of internal junctions, Xiji^)- 

The phase difference in the leads is fixed by the applied voltage bias, Xr ~ Xl = 2eVt/h, giving rise to a steady 
non-dissipative current. We will take into account the interaction with the thermal bath in the leads by adding the 
Langevin thermal force tpit) generating phase fluctuations in the leads. This means that the effective constraint on 
the phase on the leads can be written as: 

Xn-XL = 2eVt/h + ^it) , (15) 
is determined by the classical Nyquist noise in an external 



where the phase correlation function in the leads Ki^ads 
resistanceHi65^. 



Kieads = (exp['0(t) - '0(O)])r = exp 



duj 



Re 



coth 



[coB{ujt) — \\ — i snv{ix}t) 



(16) 



Thus the correlation function Kif,ads determines the interaction of the JJA with the thermal bath, and the energy 
relaxation takes place only in the leads, but not inside the array ((...)t stands for thermodynamic average). Shifting 
all internal phases Xij ~ Xm as Xij ~ Xm — ^'^ t + Xij — Xu, we bring Hq to the form (we omit the tilde-sign): 



Ha = 



[j^ ^X^ J -Xkif - E. J cos{^ 



2eK: 



ij — kl ^ 



+ X^] -Xkl) \ 



{i].kl) 



cO 



(17) 



where Vij^ki{t) are voltage drops between the adjacent islands. Plugging ([TSj) into (I13|) we find that the term H* 
in the Hamiltonian can be viewed as a time-dependent perturbation (the ac Josephson current) oscillating with the 
frequency ui = 2e(Vij -|- VNj)/fi: 



M 

Ej ^ cos 



i^it)-Xij{t)+XNAi) 



(18) 



and that the average value of the dc current can be con- 
sidered as calculated in the first order with respect to 
this perturbatioii^iiZI. Here, Vij = Vi and Vnj — 
are the dc voltage drops between the left lead and first 
row of islands, and the iV-th row of islands and the right 
lead, accordingly. Following the general recipe^^i^liZ^, we 
carry out the averaging in Eq. with the help of the 
nonequilibrium density matrix p(i), which satisfies the 
equation 

m = ~iL{t)p{t) , (19) 

where the Liouville operator is determined as 

LX = \ [Ho + Hbath + mnt+H*,X]. (20) 



Solving the Eq. up to the first order in H* we obtain 



Pit) -PP-Jil dse-^"^'-^^ [H\ Pp] , (21) 



where is the equilibrium density matrix correspond- 
ing to the Hamiltonian iJp + Hbath + Hint , the Liouville 
operator Lq is determined as LqX — {l/h)[Ho + Htath + 
Hint J X] . The expression for the average dc Cooper pair 
current assumes the form: 



, , AEj a^E"] 1 



dt I ds sin 



2eV 



(t-s) 



[Fi{t-s),Fi] + [F2{t-s),F2] 



(22) 



Ho-.i: 



where [Fi,2 {t — s), ^1,2] is the commutator of the operator Fi,2 and the corresponding Heisenberg operator F1.2 {t — s). 
The functions Fi 2 are determined as 



Fi = cos [ip - Xl] + XNj] , F2 = sin [V' - Xij + XNj] 



(23) 
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We define the time-dependent correlation function Ktot{t) of a whole system (including leads) as 

Ktot{t) = (^expS^i[m-m-Xidt)+XidO)+XNj{t)-XNj{0)]'^^^ ^ ■ (24) 

Using the property of average values^i (F(t)G'(0))p^ ~ {F{0)G{—t)) and expressing all the commutators through 
Ktot{t) we arrive at the general equation for the dc current 

Qm[Ktot{t)] . (25) 
Since Hq does not contain the variable ijj and Hint: its turn, docs not depend on the intrinsic variables Xiji 

the 

correlation function ((24|) factorizes: 

Ktotit)=Ki,ads{t)K{t), (26) 

where the correlation function of the phase noise in leads is determined by Eq. (fT6|) reflecting as we have already 
mentioned the interaction with the thermal bath (see ([5]) and ([6])) and 

Kit) = (expi[xij{t) - XI, (0) - XNjit) + XN,m)Ho , (27) 
is the time-dependent correlation function of the intrinsic p art of the system with Hp defined in (fTT]). 



dt sin 



Note here that in a two-junction system (a single 
Cooper-pair transistor) where xij = Xnj and K{t) = 1, 
the current-voltage characteristics I{V) is determined by 
the external resistance^°. We calculate the I-V curve for 
the current flowing through a single Cooper-pair transis- 
tor as an example of an application of the general formal- 
ism. If the external resistance Rext < Rcpi the current 
displays a peak in the low-voltage region (this state is of- 
ten referred to as the "pseudo-superconducting state"): 



I{V) ~ MhE,, 



Ec 



Re 



eV 



+ {kBTRext/Rcpf 



(28) 

Proportionality of I{y) to E^j indicates the Cooper-pair 
cotunneling type of transport. Note, that the Cooper- 
pair current in the two-junction system does not dis- 
play Coulomb blockade effect when the external resis- 
tance Rext < Rcp- In the opposite case, Rext > Rcpi 
the Eq. with K{t) = 1 yields the Gaussian IjV) de- 
pendence / cx exp[-(ey - E^f/^kBTEc)] (see Ref. 

Now we discuss the general case of a large size JJA 
with the number of Josephson junctions larger than two. 
We consider the case of the small external resistance, 
Rext < Rcpi which is the most frequent experimental 
situation. Expression shows that the dc current de- 
pends explicitly upon the voltage drops Vi and Vm on 
the leftmost and rightmost junctions, while the voltage 
drops Vij-ki on the internal parts of the array come in 
through the correlation function K{t). To evaluate the 
voltage distribution along the system, we notice that in 
the insulating domain, Ej <^ Ec, the correlation function 
K{t) oscillates with the high frequency determined by the 
macroscopic collective Coulomb gap Ac 3> ksT ^ Ec 
given by Eq. and that in the time interval, from which 
the main contribution is coming from, Kicads ~ 1 (see 



two next Sections), giving rise to the exponentially low 
conductance G{T) oc Goe'^"^^''^'^^ (where Go is the 
non-activated factor in the conductance characterizing 
an individual junction), i.e. 



G{T){E,,/Ec)\Vi+V, 



(29) 



Then the dc current passing through an internal junc- 
tion can be estimated as [in the first order approx- 
imation with respect to the time-dependent terms 
Ej cos[{2eVij ^kit / fi) + Xij ^ Xki] in the Hamiltonian Hq 

of uni)]: 



(int) 
dc 



GoVij-kiiEj / Ec)'^ 



(30) 



[cf Eqs. (I2ip - (P5|) ]. Utilizing the dc current conservation 

law (the Kirchhoffs law), i.e. the fact that / ~ Idc*^, we 
find: 



V 



Vi+V„ ~ ^ 

l+[NE^,/{E,EcY] [G(T)/Go] 

and, therefore, as long as 

{E.,Ecf Go 



, (31) 



(32) 



the highly inhomogeneous voltage distribution takes 
place, i.e. almost all the applied dc bias V drops on 
the first and the A^-th rows of junctions, Vi +Vn ~ V. 
In this regime, the dc current / ~ G{T){Ej/Ec)'^V 
fiowing through the array seems not to depend on the 
Josephson coupling of intrinsic Josephson junctions Ej. 
However, as Ej becomes too small and/or, talking about 
the 2D case, the number TV of junctions is growing too 
large, the condition ([5^ breaks down. This gives rise 
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to the even distribution of the total voltage drop along 
the whole array. Yet due to extremely small values of 
G{T) in the insulating regime, there exists a wide range 
of parameters where the highly inhomogeneous voltage 
distribution giving rise to the synchronized collective be- 
havior is realized. 

The above estimates suggest the following simple pic- 
ture of the macroscopic Coulomb blockade governing the 
Cooper pair insulator dynamics: the applied voltage dis- 
tributes mostly between the leftmost and rightmost rows 
of junctions, while the internal part of an array acts as 
a coherent superconducting island providing the macro- 
scopic Coulomb barrier Ac ~ (2e)^/2Ctot, where Ctot is 
a total capacitance of the array. 



III. TRANSPORT IN JJAS: THE INSULATING 
REGIME 

In this Section we calculate the current defined by Eqs. 
psp - p7p . in the insulating state, Ec ^ Ej, at moder- 
ately low temperatures Ac > ksT > Ec- 

To begin with, we find the correlation function K(t) 
for a simplest ID case. In the first approximation we put 
Ej — 0, and, to simplify the notations, assume Eco = 
oo (the generalization to the finite Eco case is almost 
straightforward). Then the Hamiltonian Hq = 
where Hi are the Hamiltonians of individual Joscphson 
junctions. Using the quantum-mechanical definition of 
K{t) we obtain^! 



K{t) 



t[Err,-Er,] 



N 



where E^s are the energy levels of a single junction. 



These energy levels are determined by the charge energy 
i?c as En = {Ef.'n?)/2, {n\e'^f\m) are the matrix elements 
of the operator e^^ between the n-th and m-th states, 
and A = exp(— E'^/fesT)]"^ is the normalization co- 
efficient (which cancel out from the final expression for 
K{t)). Therefore, the quantum mechanical dynamics of 
a single junction can be mapped onto a well studied be- 
havior of quantum rotator which has the matrix elements 
for the operator e'*^ between the states n and n + 1 only. 
We find 



Kit) 



inEct- 



N 



At Ec < ksT we replace the sum by the integral 



Kit) = e»^^=*/2 



A J dx exp{iEctx - Ecx'^ / {2kBT)} 



N 



and arrive at 



Kit) 



_ iNE^t/2~NE^kBTt^ /2 



(33) 



In a general ID case with the finite value of Eco the 
time dependence of the correlation function K{t) pre- 
serves its form ((33)) but the quantity NEc has to be 
replaced by a more general expression determined by 
the full capacitance matrix of the array. A crucial con- 
ditions allowing to obtain this result are the bilinear 
form of the Hamiltonian Hq in the momentum repre- 
sentation, and the factorization of the wave functions: 
\n)Ho = nr=i(l/N/27i^)exp(i0fenj). 

In the two-dimensional situation the calculations are 
more involved and the correlation function of the array, 
K{t), is derived as an analytical continuation of the quan- 
tity K{t), where r is the imaginary time. 



K{t) = J Li[;^y]e'[xi.(^)-xi.(o)-Xiv,(r)+xiv,(o)] 1^ _ _ y 



S e: 

{ij.kl) 



(34) 



~Ej cos[x^j (t) - Xkl (t)] 

Note that the correlation function K{t), determining the ac synchronization between the external leftmost, 1-st, and 
rightmost, iV-th, Joscphson junctions contacting with the left and right leads respectively, is not zero even in the 
zero-approximation, Ej = 0, with respect to the intrinsic Joscphson coupling Ej. The phases Xij ^-^e written as 
Xij + 2t: MijiksTr /H), where Xij is a periodic function on the interval < r < h/{kBT), and My are the winding 
numbers. 

In the high temperature regime fc^T ^ Ec we can neglect all nonzero winding numbers. Indeed, the nonzero 
winding numbers contribution to the K{t) can be estimated as: 

K^-Hr) ^ exp{(M,,-M.,)^^- E ^^[M.^-M^l^. 

{ij.ki) " 

Therefore, on the time scale r <C h/ Ec the contribution of nonzero windings numbers is small. Since the characteristic 
time T ~ h/{kBT) in the integral over time in the Eq. (j25p . we neglect all nonzero winding numbers in the limit 
Ec/ksT <ti 1. The nonzero winding numbers become important at low temperatures, ksT <^ Ec- 
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Next, we expand the periodic phases Xiji''') over the Matsubara frequencies (see also Appendix): 

uj„=2TrkBTn/h 

and change the variables in the integrals over Xiji'-^n)' 

Substituting (|55)) and ([55]) into ([M)) yields the correlation function K(t) in the following form 



K{t) = J D[x.,]expl^J2 



ErXj 



(ij.ki) 



(35) 



(36) 



(37) 



The function K{t) is the periodic function of r with the 
period = h/{kBT). Therefore, it is enough to calculate 
the sum over uin in (I37p for < r < h/{kBT). In this 
range of r 



lQEckBTsin^{uj^T/2) _ 4Ec ^E^ksT ^ 



and one finally arrives at 



K{t) ~ exp 



4A, 



(38) 



where Ac is the macroscopic Coulomb gap for the Cooper 
pair propagation defined through the functional integral 
on the lattice: 



exp(-Ac/A;BT) = J D[xij]exp 



^ 2 

{ij,ki) 



[Xt^ 



Xkl 



kBT 

1 



l[X\j - XNj) 



2E, 



cO 



.(39) 



The analytic continuation of the periodic function 
K{t) to the real time t, and the corresponding calcu- 
lation of QmK{t) is to be carried out according to the 
general recipes of the statistical physics^S. First, we find 
the quantities K{ujn) determined by the Matsubara fre- 
quencies (ujn — 2'Kn{kBT)/h): 



Jo 



(40) 



[we remind that K{t) is the periodic function of t on 
the interval {0,h/{kBT))]. The next step is the analytic 
continuation iujn uj + iS, giving rise to the retarded 
correlation function K^. Performing then the inverse 
Fourier transformation, we get: 



/oc 
exp {-iujt) 



duj 
2^ 



To carry out the analytic continuation we transform 
the integral in (|40p from the real axis to the contour in 



IZ 



i 



+iz 







FIG. 2: The contour for calculation of the integral in Eq. (|40 



the complex plain, i.e. two lines iz and h/{kBT) + iz, 
where z runs first from cx) to and then from to oo 
(see Fig. [2]) and find 



I I dz 
'o 



K{iz)~K{iz+^) 



(41) 



After that, wc change iujn to + i5, and performing the 
inverse Fourier transform, we obtain^ 



^mK{t) 



K{it)~K{it+--) 



(42) 



and in the limit Ec <C /csT <C Ac we finally arrive at 

• (43) 
we find 



c^„.X(t).exp(-1^4f^^sin^4^^* 



Calculating the integral over time in Eq. 
the I-V dependence in the following form 



(2Ac - e{Vi + y„))- 

4fcBrAc 



(44) 



(V^i + V;v)exp 



Substituting the expression for Vi + V,v ([31]) in (j44|) we 
obtain the I{y) characteristic of a large JJA in a form: 



I (X V exp 



(2Ac - eVf 
^kBTA, 



(45) 



with 



{EjE^Y Go J 



9 



Note that if we set E, = in Eq. p5)) . then the current 
through the system is zero, / = 0, as it should. On the 
other hand, if the condition ([32|) is satisfied, Vi + Vn ~ V, 
and the I-V curve assumes the simple final form 



V^exp 



(2A, 



4fcBTA, 



(46) 



The Gaussian formula ([46]) was obtained earlier for a sin- 
gle JJ incorporated in a circuit with the high resistance; 
the corresponding peak in the I-V curve was considered 
a manifestation of the "Coulomb blockade of Cooper-pair 
tunneling"—. Experimentally such a peak has been ob- 
served in Refs. [gtIIbH . In our case of large J J As one does 
not expect the similar Gaussian peak in the I-V depen- 
dence. The reason is that on approach of the bias eV 
to Ac, the conductance G{T) grows appreciably and the 
condition (|32p breaks down. The voltage distribution be- 
comes homogeneous and formula (j46p does not hold any 
more. 



At low voltages, eV < Ac, Eq. (gH]) yields the ther- 
mally activated behavior of the resistance: 



Ra 



oc exp 



Ac 



(47) 



and therefore we identify the activation temperature Tq 
of the experiment with the macroscopic Coulomb block- 
ade barrier Ac/fcs. In order to carry out calculations 
in Eq 

"spin- wave types" fluctuations: 



39|) and determine Ac, we consider the standard 



dp exp {ipRij)x{p) . 



Taking all the Gaussian integrals over x{p) in Eq. ([39]) we 
obtain the expression for Ac in the following form (the 
vector L is directed along the current in J J A): 



2Erd" 



sm 



2 pL 



(2^)" [p2 + (1/Ac)2 



(48) 



where Ac = d^J Ecol Ec is the correlation length in the 
charge coupled tunnel junction arrays, and n — 1,2 for 
the ID and 2D JJAs, correspondingly. For a large one- 
dimensional array, L 3> Ac, 



Ac = E,\J{2d) . 



(49) 



In the opposite limit of shorter one-dimensional arrays, 
i.e. L <^ Xc, the Coulomb gap increases with the array 
size L linearly: 



Ac = E,L/{2d). 



(50) 



In two-dimensional junction arrays the Coulomb gap ac- 
quires the logarithmic form: 



Ac = Ec In 



min {Ac, i} 



(51) 



IV. TRANSPORT IN JJAS: THE 
SUPERINSULATING REGIME 

Now we turn to low temperatures, fcsT <C Ec, where 
all the windings numbers Mij have to be taken into ac- 
count. In a one-dimensional array the calculation of the 
correlation function K{t) is straightforward. Consider, to 
be specific, the case where the screening length is larger 
than the size of a system L, then: 



Kit) 



N 



A = 



(52) 



and the values n = 0,±1 give the main contribution. As 
a result we arrive at 



K{t) = exp 



iAct-f2Are-^<=/(2fei5i^)[cos(£:ct)-l]| . (53) 



Substituting ([53| into Eq. ([25| we find the expression for 
the current as 



I{V) cx exp 



{eV - Ac)2e^=/2fc«^ 



8E?N 



(54) 



This result holds in the temperature range Ec/ In(iV) < 
ksT « Ec. 

In the two-dimensional case, the contribution from the 
nonzero windings numbers is analogous to the vortex con- 
tribution which appears in the classical two-dimensional 
planar Hciscnberg model (or classical Josephson junction 
arrays below the BKT). We will follow the procedure 
developed for calculation of the vortex contribution to 
various correlation functions in Ref. [8l|. Namely, it was 
shown that the coordinate dependent correlation func- 
tion 

where r and r' are the two points on the 2D lattice, can 
be expressed as 



5('„orte.)(r-r') ^ exp (-V^ In- 



where ^ — X^ro ^o('^(0)"^('"o)) is the space correlation 
function of vortex (charge)-antivortex (anticharge) pairs, 
diverging near the binding-unbinding transition temper- 
ature. Using this result and taking into account the cor- 
responding mapping p = rEc/h, we find 



This concludes the description of thermally activated be- 
havior in the temperature interval Ec < ksT < Ac. 



K{t) = exp - 



AcEc^t^ _ . 2 Act 
fi2 ' h 



(55) 
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Plugging ((55)) into (|25|) and calculating integral over time, 
we obtain, at low voltages, the following expression for 
the resistance of the array: 



Superconducting island 



cx exp 



(56) 



At low temperatures, fcsT <^ Ec the concentration 
of the charge-anticharge pairs is small and accordingly, 
C = const • exp[~Ec/{kBT)] (see Ref. jSli])- This gives 
double-exponential behavior of the resistance 



R cx exp 



A, 



exp 



(57) 



in the superinsulating regime. We would like to empha- 
size here that the double-exponential temperature depen- 
dence favors enormously the fulfilling the condition (j32[ 
for the inhomogeneous distribution of the voltage drop 
which ensures the validity of our approach. 



V. A QUALITATIVE PICTURE 

To gain physical insight in the transport phenomena 
near SIT of the large one- and two-dimensional Joseph- 
son junction arrays and films, let us discuss the distri- 
bution of the electric field in the experimental systems 
in question. Consider first one- and two-dimensional 
JJAsiiSii^ii^i^i^iiiSii^. The arrays are comprised of over- 
lapping superconducting platelets (islands) separated by 
thin oxide layers (see Fig. [3]). The related junction capac- 
itance, C, well exceeds the capacitance of each constitu- 
tive island to the ground, Cq. Thus, the total capacitance 
of the JJA is determined by the capacitance of a junc- 
tion C . If now we place a charge in the array, the induced 
electric field will remain within the array plane. In other 
words, one- and two-dimensional arrays can be viewed as 
systems with the anomalously large dielectric constant 
£ ~ C /Cq. Accordingly, in ID arrays charges interact lin- 
early over distances Z < Ac, and in 2D arrays the charges 
interact logarithmically over scales I < Ac Ref. [3]. 

Turning to disordered superconducting films and gran- 
ular superconductors near the SIT, we recall that in the 
two phase system in the vicinity of the percolation tran- 
sition between the conducting and insulating phases, the 
dielectric constant diverge c^^i^'^ . The origin of this diver- 
gence can be understood from the simple picture of the 
percolation transition, see Fig. 21 The white spots repre- 
sent superconducting clusters in the insulating sea (dark 
blue). The capacitance between the two adjacent clus- 
ters is proportional to the length of the insulating layer 
separating them. Upon approaching the transition from 
the insulating side of the SIT (see Fig.|3|D), the length 
of this layer diverges infinitely. It results in the diver- 
gent growth of the effective capacitance of the system, 
implying the divergence of the dielectric constant. Since 
the SIT in disordered films and granular superconduc- 
tors is supposed to be of the percolative natur o^^'^^i^^'^'^ . 






Superconductor 


Superconductor 


e 



FIG. 3: Distibution of the electric field in one-dimensional 
array of superconducting islands connected by two Josephson 
junctions to neighbors corresponding to experimental system 
of Ref. [|. 



Superconducting 
islands (clusters) 




Dielectric interfaces 
separating So islands 



Dielectric interface 
separating "left" and "right" 
superconducting clusters 
at the percolation threshold 



FIG. 4: Developing of the large capacitance in the disordered 
film on approach to the superconductor-to-insulator transi- 
tion from the insulating side, (a) Superconducting clusters 
separated by dielectric interfaces, (b) The same picture at 
the immediate vicinity of the SIT. Cutting the last interface 
separating the left and right superconducting clusters (green 
line) gives rise to percolation over superconducting areas from 
the left electrode to the right electrode, i.e. to the supercon- 
ducting state. 



one expects that these systems possess anomalously large 
£ near the transition. Recently the enhanced dielectric 
constant was indeed observed near the SIT in ultrathin 
amorphous beryllium films^. We therefore can conclude 
that superconducting films near the SIT exhibit two di- 
mensional behavior with respect to Coulomb interaction. 
Note that from the viewpoint of the Coulomb interac- 
tion between charges, the 2D JJAs and disordered films 
near the SIT are alike: what matters is the logarithmic 
interaction between the charges, while the differences in 
internal structure between the systems in question are 
irrelevant (at least in the absence of the magnetic field). 
The two-dimensional character of the Coulomb inter- 
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action has important implications for the systems with 
the size less then the electrostatic screening length Ac- 
Unbinding of logarithmically interacting topological exci- 
tations gives rise to the celebrated universal Berczinskii- 
Kosterletz-Thouless (BKT) phase transition first intro- 
duced in the context of vortices in XY -magnets and 
extended then to vortex-antivortex pairs in superfluid 
and superconducting films and Josephson junction ar- 

„„„ „84.85.86.87 

raya ' i i . 

On the superconducting side of the superconductor- 
insulator transition logarithmic interaction between vor- 
tices gives rise to BKT transition separating the super- 
conducting low-temperature phase, where vortices and 
antivortices are bound in pairs, from the high temper- 
ature phase with free vortices. In the high-temperature 
domain the free motion of vortices breaks down the phase 
coherence and a superconductor falls into the resistive 
state with the resistance much less than that in the nor- 
mal state. At T = Tg^T — Ej/kg the phase coher- 
ence restores and the 2D array or film becomes super- 
conducting. On the insulating side the film experiences 
charge binding-unbinding transition at the temperature 
T = Tsj~ Ec/kB [iliOIillillli, dual to the BKT in the 
superconducting state. In the high temperature phase, 
T > Tsi the charges of either sign form a gas of free 2e 
charge Cooper pairs. At low temperatures, T < Tgi, the 
charges of the opposite signs are bound in dipoles. This 
charge binding-unbinding BKT is a realization of the 
earlier theoretical observation that in a two-dimensional 
electrolite with the logarithmic interaction between ions, 
the transition at which the ions of the opposite sign be- 
come bound in pairs, occurs upon lowering down the tem- 
perature of the systemSS,. 

Having established a background, we are now in a po- 
sition to give a qualitative picture of the transport in a 
large JJA. The electric field induced by the charge placed 
on a superconducting island (or distributed over several 
islands) remains trapped within the JJA. Thus a JJA is a 
one- or two-dimensional system with respect to the elec- 
tric field distribution. Let us consider the situation where 
the screening length, Ac, appearing due to capacitance to 
the ground exceeds the sample size L. In this case the 
energy necessary to place an additional Cooper pair into 
the system is Ac = Ec{L/d) in one- and Ac = Ecln{L/d) 
in two-dimensional case. This Coulomb energy can be 
presented as Ac = (2e)^/2Ctot, where Ctot is the total 
capacitance of the array. In an one-dimensional regular 
array Ctot = C/N (the total capacitance for the system 
comprised of N capacitors in series). Correspondingly, 
in a two-dimensional system Ctot — C/\nN for the ar- 
ray containing Nx N junctions. The thermally activated 
transport is thus governed by the activation barrier Ac 
and the corresponding resistance is 



R (X exp(Ac/T) , 



(58) 



vated resistance exists only in a moderate temperatures 
region Ac > fc^T > Ec, above the transition temper- 
ature Tsi, where the free charges can propagate across 
the array. To understand the dynamics below the Tsi, 
let us notice first that the thermally activated behavior 
(j58p means that the whole array acts in a synchronized 
manner as a one single superconducting island with the 
characteristic capacitance Ctot- In the insulating state 
Ec ^ Ej the charges at every junctions are fixed and thus 
by the quantum mechanical uncertainty principle the cor- 
responding phases fluctuate loosely and so do the related 
local electric fields. However, as the current starts pass- 
ing the array, all the internal phases synchronize in order 
to minimize the Joule losses. Thus the phase evolves co- 
herently over the array implying that the whole system 
behaves as a single superconducting island. 

As a next step, one can realize that the Cooper pair 
propagation across the array can be viewed as a propa- 
gation of a charge soliton which is not necessarily con- 
fined to a one island. Following— we introduce the lo- 
cal charge density, ns(r), which is normalized to give 
the total soliton energy as Ac = Ec J drnKr). The 
probability for such a local density to appear at point 
r is proportional to exp[— ng(r)/(2(5n^))], where {Sn'^) 
is the mean square fluctuation of the local charge den- 
sity. The Cooper pair current is proportional then to 
the number of solitons generated per unit time and 
traversing the array. The latter is proportional to a 
product of all the above local probabilities at all the 
points of the system: cx Y\^exp[—nl{r)/{2{6n'^))] = 
exp{-[l/(2(5n2))]/drn2(r)} = exp[-Ac/{2Ec{Sn^))]. 
At temperatures above the charge binding-unbinding 
transition, T^j ~ Ec/kB, the solitons are unbound 
and, according to the equipartition theorem, (Sv?) = 
kBT/Ec, giving rise to thermally activated resistance 
R (X exp(Ac/(A;BT)). At low temperatures, T < Tsi, the 
charge solitons and antisolitons are bound, and there- 
fore (Sn?) is the probability of breaking these pairs, i.e. 
exp{~Ec/ {ksT)). This yields a double-exponential resis- 
tivity in the superinsulating phase: 



R cx exp 



E, 



exp 



Ec 
knT 



(59) 



reproducing the experimentally observed dependence (1) 
in the Cooper pair insulating state. The thermally acti- 



The transition from the thermally activated insulating 
to superinsulating behavior can be viewed as a manifes- 
tation of the fact that {5n?) represents the mean filling 
density n for the energy state E = Ec hy the Cooper 
pairs. The filling density n, in its turn, is given by the 
Bose statistics: Sn^ = n— [exp{Ec/kBT) — 1]~^. 

The outlined picture of the Cooper pair transport im- 
plies that the internal part of an array acts coherently 
as a single superconducting island, while the most of the 
applied voltage drops at the leftmost and rightmost junc- 
tions. In other words the system can be viewed as a 
two-junction system with the capacitance between the 
central island and leads equal to the total capacitance of 
the array. The criterion for this scenario to hold at tem- 
peratures T > Ec/kB can be presented as [see Eq. ([5^ 
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above] : 



N 



exp 



Ac 



< 1. 



(60) 



One sees that in a one-dimensional system, where Ac ~ 
NEc, the larger the system, the better the criterion ([50)1 
is satisfied in compliance with the experimental observa- 
tion o^ that the insulating behavior of the ID Josephson 
array becomes more pronounced with the increase of the 
system length. In the 2D case the situation is more com- 
plicated, but this criterion is satisfied pretty well at low 
enough temperatures ksT > E,,. In the superinsulating 
phase, ksT < Ec, the corresponding criterion following 
from Ea. (j32p is met very well. 

We expect that the model of the large Josephson jimc- 
tion array applies fairly well to thin films in the criti- 
cal region of the superconductor-to-insulator transition, 
since as we have already noticed, the internal structure 
is irrelevant with respect to Coulomb properties of the 
system. 

In conclusion, we have shown that the Cooper pair 
transport in the insulating state of one- and two- 
dimensional Josephson junction arrays is governed by 
the macroscopic Coulomb blockade. The macroscopic 
Coulomb blockade energy Ac ~ Ec{L/d) in ID and 
Ac ~ Ec\n{L/d) in 2D systems. We have shown that the 
charge binding-unbinding BKT-like transition separates 
the insulating high temperature state with the thermally 
activated conductivity from the low temperature superin- 
sulating state. We have determined the conditions under 
which the macroscopic Coulomb blockade is realized and 
the conditions under which the Coulomb blockade acti- 
vation energy exhibits the system size dependence. 

The questions that remain open include: 

1. The contribution of the quasiparticle current into 
the transport properties of large Josephson junc- 
tion arrays. 

2. The role of quantum fluctuations in the transport 
properties of large Josephson junction arrays. 

These and related topics will be a subject of forthcoming 
publication. 
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APPENDIX 

As an illustration, we calculate the correlation function 
for a single JJ case with the Hamiltonian 



where ip is the Josephson phase. In this case, 

l-h/(kBT) 



K{t) = J D(^e*['^(^)-'^(o)l exp 



4-Bc 



drip 



The calculation is done via expanding ip{t) into a Fourier 
series: 



ip{t) 



E 

^=2TrkBTn/h 



Replacing the functional integration over ip{t) by inte- 
gration over Fourier coefficients ipn yields 



fiS, ,2, .2 



^E^keT 



-|-^„(e-"--l) 



^EcksT . 



In the interval < r < h/ [ksT) the sum over n yields: 



E 



4£;cfcBrsin'(^nT/2) _ E, E.keT _ 

and the correlation function assumes the form 
if(r)^exp(^- — + 

Finally, the changing r to it gives 

( .Ect E.ksTt^ 

This result coincides with the one obtained by direct 
calculations using the quantum-mechanical definition of 
K{t) (Ref. [Zl) and 
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